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Abstract 



The problem of estimating a complex measure made up by a linear combination 
of Dirac distributions centered on points of the complex plane from a finite num- 
ber of its complex moments affected by additive i.i.d. Gaussian noise is considered. 
A random measure is defined whose expectation approximates the unknown mea- 
sure under suitable conditions. An estimator of the approximating measure is then 
proposed as well as a new discrete transform of the noisy moments that allows to 
compute an estimate of the unknown measure. A small simulation study is also 
performed to experimentally check the goodness of the approximations. 

Key words and phrases: Complex moments; Pade' approximants; logarithmic po- 
tentials; random determinants; random polynomials; pencils of matrices 
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Introduction 



Let us consider the complex measure defined on a compact set D C (T by 

S{z) = T.c,5{z-Q, eint(D), c.eW 

and let be 

= J z^'S{z)dz = JJ{^ ^ iy)^S{x + iy)dxdy, A; = 0, 1, 2, . . . 

D D 

the complex moments. It turns out that 

sk = j: c.e,. (1) 

Let us assume to know an even number n > 2^ of noisy complex moments 

^k = Sk^h'k, A; = 0, 1, 2, . . . ,n - 1 

where Uk is a complex Gaussian, zero mean, white noise, with finite known vari- 
ance cr^. In the following all random quantities are denoted by bold characters. We 
want to estimate S{z) from {afc}fc=o,...,n-i- From equation ([T]) this is equivalent to 
estimate j9, Cj, ^j, j = 1, . . . , j9, which is the well known difficult problem of complex 
exponentials approximation. 

The problem is central in many disciplines and appears in the literature in different 
forms and contexts (see e.g. [IH lfT^lE^lEHES] ) . The assumptions about the noise vari- 



ance (constant and known) are made here to simplify the analysis. However in many 
applications the noise is an instrumental one which is well represented by a white 
noise, zero mean, Gaussian process whose variance is known or easy to estimate. A 
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typical example is provided by NMR spectroscopy (see e.g. jB]). 

In the noiseless case the problem becomes the complex exponential interpolation 
problem [H]. Conditions for existence and unicity of the solution are ([HI Th.7.2c]): 



detUo{s) 7^ 0, detUi{s) 



where 



U{so, . . . , S2p-2) = 



So Si ... Sp^i 



Si S2 



Sp—1 Sp 



S2p-2 



and 



^o(s) = U (so, . . . , S2p_2), Ui{s) = U (Si, . . . , S2p_l). 

In fact exactly n = 2p noiseless moments are sufficient to fully retrieve S{z), where 

p = max{n | det{U{so, . . . , s„_2)) 0}. 

Moreover (^j, j = are the generalized eigenvalues of the pencil P = 

[Ui{s), Uo{s)] i.e. they are the roots of the polynomial in the variable z 

det[Ui{s) - zUo{s)] 

and Cj are related to the generalized eigenvector Uj of P by cj = uJ[so, . . . , Sp_i]^. 
In fact from equation ([T]) we have c = V^^fso, . . . , Sp-_i]^ where 



V = Vander{^i, . . . ,^p 



is the square Vandermonde matrix based on (^i, . . . But it easy to show (see 
e.g. [2j) that 

Uo{s) = VCV^, Ui{s) = VCZV^ 

where 

C = diag{ci, . . . , Cp} and Z = diag{^i, . . . , ^p}. 

Therefore Uk = V^^Ck is the right generalized eigenvector of P corresponding to 
where is the k—ih column of the identity matrix /„ of order p. 



Viceversa when Sk = 0, V/c it was proved in |T5] that 

det[U{suQ^ . . . , a„_2)] = det[UQ{a^] 7^ Vn a.s. 

and 

det\U{aii, . . . ,a„_i)] = det\Ui{^] 7^ Vn a.s.. 
Moreover associated to the random polynomial 

det[Ui{^) - zUq{^)] (2) 

a condensed density hn{z) can be considered which is the expected value of the 
(random) normalized counting measure on the zeros of this polynomial i.e. 

2 



hn\z) = -E 
n 



n/2 

E Kz - 



It was proved in [1] that if 2; = re*^, the marginal condensed density h^^\r) w.r. to 
r of the generalized eigenvalues is asymptotically in n a Dirac 5 supported on the 
unit circle Vcr^. Moreover for finite n the the marginal condensed density w.r. to 6 is 
uniformly distributed on [— 7r,7r]. Starting from the generalized eigenvalues $,j and 
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generalized eigenvectors Uj of the pencil 

P = [t/(ai, . . . ,a„_i), t/(ao, . . • , a„_2)] 
we then define a family of random measures 

n/2 
J = l 

where Cj = uJ[ao, . . . , a„/2-i]"^ and we give conditions under which E[Sn{z)] approx- 
imates S{z). Moreover we define a discrete transform (P-Transform) on a lattice of 
points on D, which is an unbiased and consistent estimator of E[Sn{z)] on the lattice 
thus providing a computational device to solve the original problem. 

In jl] the same problem was afforded. The joint distribution of the coefficients 
of the random polynomial (Ej) (when Sk 7^ 0, V/c) was approximated by a multi- 
variate Gaussian distribution and a theorem by Hammersley [7j was used to com- 
pute the associated condensed density of its roots. An heuristic algorithm was then 
used to identify the main peaks of the condensed density and to get estimates of 
j9, and Cj, j = 1, . . . ,j9 based on them. In the present work the ideas presented in 
[Ij are put on a more rigorous mathematical framework. A different approximation 
of the condensed density is considered and an automatic estimation procedure is 
proposed. 

The paper is organized as follows. In the first section we study the distribution of the 
generalized eigenvalues of the random pencil P and we give an easily computable 
approximate expression of the associated condensed density. In section 2 we consider 
the identifiability problem for S{z) given the data a. Conditions for identifiability 
are given and the approximation properties of E[Sn{z)] are proved. In section 3 the 



6 



P-transform is defined and its statistical properties are studied. In section 4 the 
procedure for estimating the parameters j9, {^j,Cj,j = of the unknown 

measure from the P-transform is described. FinaUy in section 5 some experimental 
results on synthetic data are reported. 



1 Distribution of the generalized eigenvalues of the pencil P 



We start by making some technical assumptions on the noise model. When Sk = 
V/c, we noticed in the introduction that are, asymptotically on n, uniformly 
distributed on the unit circle. Therefore, when 7^ is given by (II]), we can assume 
that Up = n/2—p among the ^j, j = 1, . . . , n/2 are related to noise and then they can 
be modeled for large n by = e "p i.e. by uniformly spaced deterministic generalized 
eigenvalues. Therefore the Vandermonde matrix based on ^j^j = 1, . . . , is simply 
given by V = y/n^ ■ F G (T^p^^p where Fhk = -h 
transform matrix. Hence 



27Tihk 



is the discrete Fourier 



1 



Hi 



V % 

and c has a complex multivariate Gaussian distribution with 

a' 



E[cj] = and E[cjCh\ = — Sjh- 



rir 



Based on these observations we define a new noise process as 



t^k = s 



^k, 



k > fir 
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and we assume that c is independent of i/^, /c > Up. But then E[i>]^] = and 



-1 ^ ~/i o ZTVvrytz— n ) 

E-f iHjE[ciCj] = ^ jfrLi = k,h< 



27rir[k—h) 



rir 



h > Up^k < Up 



h^k > rir 



We have then proved the fohowing 

Lemma 1 The random vectors i/^ and D^, k = 0, . . . ,n — 1 are equal in distribu- 
tion. 



As a consequence in the following we will use i^k without loss of generality. 



Remark 1 We notice that when Sk 0, if the signal-to-noise ratio is defined as 
SNR = ^mmh=i,p \ch\ we have 

PUp. |2l _ _ ^m=l,p \Ch? 

^' ^ rip ripSNR^ ' 
If SNR > then E[\cj\^] < \ck\^ ^j,k. 

A basic result which will be used extensively in the following is given by 

Lemma 2 Let T = {T^^\T^'^^) be the transformation that maps every realization 
a(0) of a to (^(0), c(0)) given by 0^(0) = ^]=i Cj(0)^j(0)^, A; = 0, . . . , n — 1, where 
0^0, and Q, is the space of events. Then T is a.s. one-to-one. Moreover, for cr — ^ 
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and for j = 1 , . . . , n/2 



ij-p + o{a), j =p+l,...,n/2 



E[c,] = 



o{a), 



j =p+l,...,n/2 



proof 



From [15] we know that a.s. detlUhil^)] 7^ 0, /i = 0, 1. Moreover, with probabihty 
1, there is no functional dependence between and s. Therefore a.s. det[Uh{s,)] ^ 
0, /i = 0, 1. But then a.s. the complex exponential interpolation problem for a has 
an unique solution V0 hence T is a.s. one-to-one. The second part of the thesis is 
based on a Taylor expansion of T around a suitable point xq. A natural candidate 
for Xq would be s. However we notice that T^^\s) is not defined if n > 2j9, and, as 
a consequence, also T^'^\s) is not defined in this case. Therefore, by using Lemma 
[TJ without loss of generality, we assume that the noise is represented by Ok i.e. 



afc = < 



n 



where Up = n/2 — p. We then define a new sequence Sk by 

p n/2 

j=i j=p+i 
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and we consider the process as a perturbation of Sfc. Therefore we choose xq 
and notice that 



= s 



r(2)(l), = 



c^", j = 13+ 1, • • • ,^/2 
We now prove that each component of T'^^\a) is an analytic function of a when a 
belong to small neighbor of s. The proof follows closely [27j [Th.6.9.8]. For each fixed 
0, the polynomial 

4>{z,a) = det[Ui{a) — zUo{a)] 
is an analytic function of z and a. Let ^ be a zero of 0(2;, s) and 



K={C\\C-i\=r}, r>0 

be a circle around ^ not containing any other generalized eigenvalue of the pencil 

P = [U{si, . . . ,s,j_i), t/(so, . . . ,S„_2)]- 

We want to show that K does not pass through any zero of 0(0,a). In fact by the 
definition of K it follows that 

inf |0(C,s)|>O. 

But (l){z,a) depends continuously on o, hence there exists B = {x G (r^\\x — s] < 
p}, p > such that 

inf |0(C,a)| > 0, Va G B. 



10 



By the principle of argument, the number of zeros of (f){z,a) within K is given by 

,w X 1 / (l)'{z,a) d(j) 

which is continuous in B] hence 

1 = N{1) = N{a), aeB. 
Moreover the simple zero ^(0) of (J){z^a) inside K admits the representation (see e.g. 

m) 

27rz ^ (p{z,a) 

For a E B the integrand is an analytic function of a and therefore also ^(0) is an 
analytic function of a when a E B. 

We now consider T^'^\a). We notice that each component can be obtained as a 
rational function of the components of T^^\a) by the formula Cj = ejy^^a, j = 
1, . . . , n/2 where V is the Vandermonde matrix based on T'^^\a). Therefore also cj 
is an analytic function of a when aEB. 

As T'^'*) = T^^^ -^iTj'^^ is analytic for aEB, T^^^ and t)^^ are real analytic functions 
of QlRi Qli where a = + m/, (e.g. [13j [pg.99]). Therefore they admit a Taylor series 
expansion around s when aEB: 

1 Rk [a) - ^ Rk [§-) + Fffi - SffiJ + 

-i5rS(^ 

i=0 OCLli |a=l 

^ i=o j=o oaRiOaRj <~ 
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\a=s 



and analogously for Tj^'*(a). Taking expectations we get 



n/2 

E [(aiji - SRi)] = [sRi - SRj\ = cr" • d, Ci= J2 ij-p 

j=p+i 



E [{a.Ri - SRi){aRj - SRj)] = E [{a.Ri - SRi + a''Ci){aRj - SRj + a'^Cj)] 

2 



(7 

2 



and analogously for the other terms. Remembering the independence of the real and 
imaginary parts of a^, we finally get 

E[rf (a)] = rf + □ 

We start now the study of the distribution in (T of the generalized eigenvalues of 
P by making some qualitative statements already present in the literature. For 
each realization 0, let {cj(0), ^j(0)}, j = 1, . . . , n/2 be the solution of the complex 
exponential interpolation problem for the data 0^(0), /c = 0, . . . , n — 1. It is well 
known that we can then define the Fade' approximant 

[n/2 - l,n/2](z,0) = z E = Qn/2-i{z-')/ Pn/2{z-') 

to the Z— transform of {0^(0)} given by 

•00 

f{z, 0) = E ak{0)z~'' = fs{z) + fu{z, 0) 

k=0 
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where 



P oo , P C- 



fc=o j=i fc=o j=i ^ sj 

and, because of Lemma [1], 

^ c,(0) 



J = l Z - ^3 

f{z^ 0) is then defined outside the unit circle and can be extended to D by analytic 
continuation. We get then 

^n;^r'(^-<5,(0)) 



/(z,0) ftJ zq„/2-l(z)/p„/2(z) = 



and 

n/2—l p Up 

g{z,0) = \og{z-'f{z,0))= E log(2-'^,(0))-Elog(z-O)-Elog(z-O)- 

j=i j=i j=i 

We want to study the location in (J* of ^^(0). To this aim, following |T9j, we remember 
that Pn{z) = z^Pn{z^^) satisfy the following orthogonality relation 

J z^^f{z, 0)pn{z)z^dz = 0, /c = 0, . . . , n - 1 
r 

where F is a union of closed curves enclosing the poles of /(z, 0) i.e. the numbers 
^j, j = 1 . . . , j9 and ^j, i = 1, . . . , Tip. By using the Szego integral representation 
of such polynomials and a saddle point argument, it turns out that the Fade' poles 
^j(0), j = 1, . . . , n/2 , asymptotically on n, satisfy the following system of algebraic 
equations 



l,n/2 1 

or 
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l,n/2 1 n/2-1 -i 

2 E 77-7-^^^7-7+ E 



i^fc (6(0)-^j(0)) iTi (6(0) -^j(0)) 
pi 1 
-V— — r-E ^ = 0, /c = l,...,n/2 

iil(6(0)-O) hiU0)-Q 

These equations can be interpreted as conditions of electrostatic equilibrium of a 
set of charges in the presence of an electric external field corresponding to g'{z,0). 
Therefore the Pade' poles 6(0) are attracted by ^j, j = 1, . . . ,p and ^j, j = 1, . . .Up 
and they are repelled by each other and by the zeros 6j{0) of qn/2-i{z)- However 

p l.p Up 

q.n/2~i{z) = E q n (z - n (^ - ^k) (3) 

rip p l,np 

+ Eq(0)n(^-a)n(^-6)- (4) 

j=l k=l k^j 

As V0, |cj(0) p <C min/i | c/ip if the SNR is sufHciently high (see Remark after Lemma 
[I]), we can approximate qn/2-i{z) by 

n(^-6)E9-n(^-a) 

k=l j=l k^j 

hence Up zeros are close to ^k, and the other — 1 are close to the zeros of the 
polynomial 

p i,p 

(ip-i{z) = E Cj n(^ -^fc) 

which is the numerator of z~^fs{z). We notice that if \ch\ <C \ck\, \/k ^ h then 

i,p i,p p i,p 

qp~i{z) - E n (^ - ^k) = {z- ^h) E Cj n (^ - ^k) 

j^h k^j j=l k^j,h 

Hence, because of the continuous dependence of the roots from the coefficient of 
a polynomial, qp^i(z) has a zero as close to as \ch\ is small with respect to 
|cfc|, k h. Therefore the Pade' poles ^k{0) 
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• are attracted by , j = 1 , . . . , p 

• are attracted by , j = 1 , . . . np 

• are repelled from <^j(0), j k 

• are repelled from ^j, j = 1, . . . 

• are repelled from other p — 1 points in the complex plane which are as close to 
as \cj\ is small with respect to |c/i|, h ^ j. 

Summing up a with a large \ck\ will attract a Pade' pole without being disturbed 
by the repulsion exerted by the zeros of qn/2-i{z)- Moreover close to such a point a 
gap of Pade' poles can be expected because of the repulsion exerted by Pade' poles 
to each other. A with a small \ck\ will still attract a Pade' pole but not so close 
because of the repulsion exerted by a close zero. The Pade' poles not related to the 
signal are expected to be attracted by which at the same time will repel them. 
Moreover they are repelled by hence they are likely to be located in between 
and far from ^k- A picture of this behavior is given in fig.l. We notice that the 
qualitative results discussed above are consistent with those obtained in [[3| under a 
more stringent hypothesis about the noise. 

We now wish to define a mathematical tool to quantify these qualitative statements. 
To this aim we remember that ^^,A; = l,...,n/2 are the generalized eigenvalues of 
the pencil P and therefore they satisfy the equation 

Pn/2(^"') = det[Ui{^) - z[/o(a)] = 0. 
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Then a condensed density hn{z) can be considered which is the expected value of 
the (random) normahzed counting measure on the zeros of this polynomial i.e. 

2 



hn{z) = -E 
n 



n/2 



The following theorem holds whose proof is the same of that of Theorem 1 in |[T|: 

Theorem 1 The condensed density of the zeros of the random polynomial Q,{z) = 
Pn/2(^"^) «5 given by 

hn{z) = -^Auniz) (5) 
where A denotes the Laplacian operator with respect to x^y if z = x + iy and 

ur,{z) = {log{\Q{z)\')} (6) 

The condensed density provides the required quantitative information about the 
distribution of the Pade' poles in the complex plane. If the SNR is sufficiently high, 
after the qualitative statements made above about the location of the Pade' poles, a 
peak of hn{z) can be expected in a neighborhood of each of the complex exponentials 
^k^k = and the volume under the peak gives the probability of finding a 

Pade' pole in that neighborhood. This is confirmed by the following 

Theorem 2 lfa>0, the condensed density hn{z,a) is a continuous function of z 
given by 



' ■^~-'^(z:'"/2-i(r"/2 
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where C = {C/i, h j} and 



7 



ifn = 2 



Moreover hn{z, a) converges weakly to the positive measure | Ej=i S{z—^j) when a 
0. 



proof 



Let us consider the transformation Tn : (C, 7) given by 

n/2 



or 



(T«{a)), = 0. (Tf (fl)).=7,- 
In the following, to simplify notations, {Tli^\a))j will be denoted by Cj{a). We have 



hn{z,a) = -E 
n 



n/2 

2 "/^ -1 

As the complex Jacobian of T^^ is (see PHTj) (n was assumed even) 



(8) 
(9) 



^c(C,7) 



7 if n = 2 

(-ir/^n;£i7,n,</.(o-a)' ifn>4 



by making a change of variables we have 
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hn{z,a) = 



niira 
2 



n/2 _^ ,2 



n/2 



n(7rcr^)"' ^ 



j-l(J<n/2-lg<n/2 



where C* = {Oi, ^ 7^ j} and 



7 if n = 2 



The integral above converges uniformly for z ^ hence hn{z) is continuous in D. 
We prove now that h2p{z, a) converges weakly to ^ Ej^i S{z — ^j) when a ^ 0. Let 
G C°° be a bounded test function supported onW. We have 

J h2p{z,a)^{z)dz 



(F 



1 p ^ p~Efc=oMij.P 
--^ ^(Q(w + i)) dy. 

As $(^) is continuous and bounded and Q is analytic in a neighbor of s by Lemma 
El by the dominated convergence theorem we get 



a->0 



1 P 



Ijm I h2p{z,a)<^{z)dz = -J2 / lim + s)) 
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TT 



2p 



1? j^i p j=i 
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because {T2l\s))j = ^j. 

Let us consider now the case n > 2p. We cannot use the same argument used for 
the case n = 2p because Cj{s) is not defined for j = ^ + 1, . . . , n/2 (see Lemma [2]). 
However by Lemma [1] without loss of generality, we can assume that the noise is 
represented by i>k i-e. 



afc = < 



Jfj^i Cj^^j +Vk: k = Up, . . . ,n - I 



where Up = n/2 — p. We then define a new process by 

p 
1 



afc = E Cj^j + ^fc, k = 0,...,n-l 



where 

n/2 

and we consider the process as a perturbation of the process a/^. Let us consider 
the pencils 

P = [t/(ai, . . . ,a„„i), [/(ao, . . • , a„_2)] 

and 

P = [t/(ai, . . . ,a^_i), t/(ao, . . • ,a„_2)]- 

We can write 

P = P + (jE 

where 

E = i[t/(0, . . . , 0, i^n,+i - rj^^+i, ^n-i - 



t/(0, . . . , 0, - r?^^, . . . , i/„_2 - r/„_2)] = [El, Eg]. 
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From IIB], in the limit for cr ^ 0, a generalized eigenvalue of P can be expressed 
as a function of a generalized eigenvalue of P and corresponding left and right 
generalized eigenvectors Vj^uj by 



^j = ^j + ^ 



vfUoUj 



= 6 + a 



e y-i(Ei-^,Eo))y-^e 



where Ug = ?7(ai, . . . , a„_i) and, by construction, 

j = '^^---^P 



= 



Cj-p, j =i?+l,...,n/2 



V = Vander{ii, . . . , 4/2), C = diag {ci, . . . , c^/2) 



and 



Vj = Uj = V 



We notice that we can write 



n/2+p 



h=i 

where 7^/1 are constants and Y/i( are i.i.d. zero mean, complex Gaussian variables 
with unit variance identified with -^[j^h — Vh]^ ^ = ''^pj ■ ■ ■ ,n — 1. 

We have 
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hn{z,a) = -E 



n 



n/2 



n 



E 5{z - 1,) 



n 



n/2 

E K^-ij) 

j=p+i 



By the same argument used for the case n = 2p it follows that h^\z^ a) converges 
weakly to | Ej=i S{z — ^j) when a ^ 0. We then consider h''^\z^ a). We have 



n 
n 



n/2 

E K^-ij) 

' ri/2 / 

E ^z-Cj.p 
j=p+i V 



— (7 



Z^/i=l Ijh Up+h 



Cj-p 



By identifying ^Cj_p, j = p + 1, . . . , n/2 with Y/j, /i = 1, . . . , n^, which are i.i.d. 
zero mean, complex Gaussian variables with unit variance, we get 



n/2 
n/2 

= E 



/ 



n/2+p 

- E IjhVnp+h - 0{a^) 



V 



(5 



(T 



Vj-P h=l 



E IjhVnp+h - 0{a ) 



(^~\V3-p\ 



dy 



TT 



n-1 



-dy', {y} = {y} \ {yj-p} 



by making the change of variable 



/TT' "-Z^+P 
w; = |j_p + ^^ E IjhVnp+h 

yj-p h=i 



we get 



TT 



(10) 



n/2+p 



/ ^ z - - — IjhVnp+h -0{a )\ dyj-p 
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Inserting this expression in ( JTUj ) we get 



/fir 



n/2 

/i(f)(^,a) = - E 

j=p+i {z - 0(a2) - <t^-_p)2 

En/2+p 
h=l ^jhVnp+h 



n/2+p ^ 



r=l 



TT' 



3-P 



dy 



and therefore 



/rir 



n/2 

J2 - 

j=p+l {z — ^j-p) 

J- f 

^^""j^n J y^p+r^ 



'3-P 



En 
k=l. 



\yk\ 



dy' = 



because 



n—l 



TT 



TT 



En/2+p 
h=l_ 'rjhVnp + h 



~J2k^i.,k^j-'P \yk\ 



dy' 



J Uup+rG -"^-dy' = 0, for a suitable hermitian matrix ^4, Vr. □ 



Remark . When the SNR is large the exponential part dominates the integrand as 
the Jacobian does not depend on a. Moreover the exponential part has relative 
maxima close to as expected. In general the integral (^!) does not admit a closed 
form expression. However when n = 2, remembering that the Jacobian with respect 
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to the real and imaginary part of a complex variable is Jr = | JcP, the integral ([Tj) 
becomes 



1 r |7-soP+l7z-nP , 1 /■ I i9 _lXz£(lii+l2£z£li 



2 



0-^(1 + bp) + bsi + SoP - If/o-^il' 

— 1^ ! LjI ! cr^(l- 

TTCT 



2q + U|2)3 



We notice that hmcr_>o ^2(^, cr) = ^(^ — si/so) = <5(2; — ^1). Moreover, when sq = 
si = we have h2{z, a) = ^^_^^^|2)2 which is independent of a^, confirming the result 
obtained in pj] for the pure noise case. 

The condensed density has an important role in the following. Therefore we look 
for an easily computable approximation. The following theorem provides a basis for 
building such an approximation : 



Theorem 3 Let be F{z,z) = (t/i(a) - zUo{a)){Ui{a) - zUo{a)) then 



E[log{det{F{z,z)})] - log{det{E[F{z,z)]}) = o{a) 



for (7^0, independently of z. Moreover 



2 

no 



E\F{z,z)\ = (t/i(s) - ^C/o(s))([/i(s) - zU^{s)) + -^A^z.z) (11) 
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where 



A{z,z] 



l + UP -z ... 



-z 1 + Ur -z 







-z l + \z\ 



proof 

let us denote by Xj the eigenvalues of F(^, z) and by /ij those of £'[F(2;, z)] , dropping 
for simplicity the dependence on z^z. Note that fij ^ E[Xj], see e.g. [5, Theorem 
8.5]. We have 

E[logidet{Fiz,z)})] = j:E[log{X,)] 

j 

and 

\og{det{E[F{z,z)]}) = Y.log{fij) 

j 

hence it is sufficient to study the difference 

E[log(A,)]-log(/i,). 

We then denote by f the vector obtained by stacking the real and imaginary parts 
of the elements {Fhk, h,k = 1, . . . , n/2) of F and consider the function 

^(f ) = log(A,) 

and its Taylor expansion around E[f\: 
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^(f) = ^(E[f]) + E^kf] {lh-E[i,]) 

h OLh 



1 



2 hk dlhdlk 
which can be rewritten as 



E[n{lh-E[iu]){lk-E[U])^ 



1 



log(A,) - log(^,) = Y.(^h{lh - E[iu]) + -Y.lhk{lh - E[h]){U - E[ik]) + 

h ^ 



hk 



and, taking expectations, 



1 



^[log(A,)] - log(/i,) = -ElhkEiih - E[h]){ik - E[l,])] + 

^ hk 



But 



F{z,z) = {Ui{s) 

-ms) 



zUois)iUiis) - zUois)) 



- zUo{m .)){Ui{u) - zUoju )) 
zUo{s){Ui{u) - zUoiiy)) 



zUo{iy)){Ui{s) - zUo{s)) 



and 



E[F{z,z)] = {Uiis) - zUois)iUiis) - zUois)) 



^E[{Ui{u) - zUo{m.)){Ui{u) - zUoiM.))] 



na 



= {Ui{s) - zUo{s){Ui{s) - zUo{s)) + ^^(z, ^) 

by a straightforward computation similar to that given in [1, Th.3] for the pure 
noise case. Therefore 



F{z,z) - E[F{z,z)] = {Ui{i^) - zUq{u )){U,{u) - zU^{u) ) 

- {Ui{s) - zU^{s){Ui{u) - zUoii/)) 



{Ui{u) - zUq{u)){Ui{s) - zUo{s)) 



na 



■A{z,z) 
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hence -^[(f/j — E[Zfi]){i]^ — E[ij^])] is a linear combination of functions of z and z with 
coefficients equal to either or a'^ because the odd moments of a Gaussian are 
zero. By a similar argument all the dropped terms in the Taylor expansion above 
will depend on even powers of a. Hence 

^[log(A,)] - log(A.,) = oia) 

independently of z,^. □ 

By noticing that |Q(2;)p = det{F{z,z)}, an approximation of the condensed density 
is then given by 

hn{z,a) = -^A logte(^)) 

where lJ-j{z) are the eigenvalues of E[F{z,z)]. Unfortunately hniz, a) is not a prob- 
ability density as it can eventually assume negative values. However the following 
results hold 

Theorem 4 The function hn{z,a) is continuous in a and in z. In the limit cases 
(7 = and {ck = 0, k = 1, . . . ,p} it is given respectively by 

^3=1 

and by 

hniz, a) = -^Awniz) 
47r 

where 

1 



Wniz) =-l0gJ2 \z\ 



2i_ 

Moreover, in this second case, \imn^oohn{z,o-) = 6{\z\ — 1). 
proof 
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hn{z,a) is continuous in a and in z because of the continuous dependence of the 
eigenvalues on the elements of the corresponding matrix. When cr = 0, let F G (T"/^'^ 
be the Vandermonde matrix such that Uq{s) = VCV^ and Ui{s) = VCZV^. Let 
V = QR be the QR decomposition of V. Then 



E[F{z,z)] = QRC{Z - zI)R^Q^QR{Z - zI)CR''Q 



But R = 



( \ 

R 



, therefore R^R = R^R: moreover Q^Q = /, hence the eigenvalues 



of E[F{z, z)] are the same of those of the matrix 



RC{Z - zI)R^R{Z - zI)CR 



H 



\ 



RC{Z - zI)R^R{Z - zI)CR'' 



H 











The non-zero eigenvalues of E\F{z^ z)] are then the same of those of the matrix 



RC{Z - zI)R^R{Z - zI)CR 



H 



We then have 



Kiz,0) = ^A E log(/i,(z)) 
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^ A log In \z - ^jf ■ \detiR)\' fl c'^ 



27rn 



op op 

— EAiog|z-^,f = -E^(^-o; 



4:7171 



n ,=1 



because logd^p) = S{z) (see e.g. [251 Pg-47]). When {c^ = 0,k = 1, . . . ,p} 
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~ 1 1 

hn{z,a) = ^Mog{det{A{z,z)}) = :^Alog(E l^f'] 



2Tin 



27rn 



The last part of the thesis follows by the same argument used in the proof of 
Theorem 3 in |jl|. □ 

Corollary 2 hn{z,a) — hn{z^a) converges weakly to when cr ^ 
proof 

Let ^{z) be a nonnegative test function supported on W. Denoting by hl^{z) = 
- E?=i 5{z — ^j), from Theorems [2] and [H we have \/v > 0, 3ai and (72 > such that 



$(2) {hn{z,a)-hl{z))dz 



< -, Vcr < (Ji 



and 



$(2) (/i„(2,(j) 



V 



< -, Vcr < (72 



hence, if Oy = minjcri, (J2}, we have Vcr < 



J $(z) cr) - hn{z, a)) dz 

(U 



< 



{hn{z,a)-h%z))dz 



w 



+ 



^{z) (hn{z,a)-h%z))dz 



<u. □ 



2 Identifiability of S'(2) and approximation properties of £'[S„(z)] 



We want now to exploit the information about the location in the complex plane of 
the Pade' poles, provided by the condensed density hn{z), to prove some properties 
relating S„(z) = J2j=iCj6{z — ^j) to the true measure S{z). 
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Before affording the problem of estimating S{z) from the data a we need to check 
that the data provide enough information to solve it. Precise conditions that must 
be met to solve the problem are well known in the noiseless case and are reported 
in the introduction. When noise is present the identifiability problem is an open 
one. Its solvability can depend on the amount of "a priori" information available 
[6] and/or on the ability to devise smart algorithms. In the following a definition of 
identifiability is given and, based on it, some properties of S„(z) are proved. 

Definition 1 The measure S{z) is identifiable from the data ak, /c = 0, . . . , n — 1 if 
3 rfc > 0, /c = 1, . . . ,]3 such that 

• hn{z) is unimodal in = {z\ \z — ^ i"k} 

The idea is that S{z) can be identified from the data a if the random general- 
ized eigenvalues have a condensed density with separate peaks centered on = 
As, by Theorem [2], hn{z,a) converges weakly to ^J2^=i6{z — ^j) when 
(7 ^ 0, it must exists a cr' > small enough to make S{z) identifiable M a < a' . 

In order to apply the proposed method one should check that the identifiability 
conditions are verified. As hn{z^ a) depends on the unknown quantities p, Cj, this is 
of course impossible. However in most real problems we have some prior information 
about the unknown measure S{z) that we can exploit to get reasonable interval 
estimates for p^cj^^j. Moreover in many instances either n or cr or both can be 
freely chosen. By Theorem El equation [HJ n should not be as large as possible to 



get the best estimates of S{z). In fact too many data will convey too much noise 
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which could mask the signal s^. We can therefore properly design an experiment by 
computing hn{z,a) for many values of n and a and choose Uou and aou (optimal 
design) that make identifiable the measures corresponding to prior estimates of 
Cj, ^j. To identify the unknown measure S{z) we then hopefully need to measure 
Uott data affected by an error with s.d. aou- Unfortunately hn{z) does not admit a 
closed form expression and to compute the expectation that appears in its definition 
we need to perform a time consuming MonteCarlo experiment. This is why we 
proposed an approximation hn{z) of hn{z) which can be quickly computed by solving 
hermitian eigenvalues problems. 

Let us consider the function 

n/2 

Sn{z) = E[Sr^{z)] = E E[CAZ - ^j] 

i=i 

where {cj, ^j}, j = 1, . . . , n/2} are the solution of the complex exponential inter- 
polation problem for the data {ak, k = 0, . . . ,n — 1}. 

The relation between Sniz) and the unknown measure S{z) is given by the following 

Theorem 5 // S{z) is identifiable from a then 

j Sn{z)dz = Ch + o(cr) 

Nh 

and 

j Sn{z)dz = o{a), V^cL»-U^j 

A j 

proof 
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From the identifiability hypothesis we know that 

I K{z)dz = - E Proh[i^ G iVfc] > 0, k = l,...,p. 

Therefore there exist such that Proh[$,j^ G Nk] > 0. Among the let us 
denote by the one such that Prob[^j^ G Nk] is maximum. From the identifiabihty 
hypothesis the are distinct. Moreover aU the j = 1, . . . ,n/2 can be sorted in 
such a way that = j = I, ■ ■ ■ ,p and, by Lemma El to it corresponds 
such that 



E[ck] 



But then for /c = 1. 



Ck + o(cr), k = l,...,p 

o{a), k = p + 1, . . . ,n/2 

n/2 

J Sn{z)dz =iZ j E[cj5{z - $,^]dz = 



'P 



n/2 



dz = 



r I r \ 



where ii^q is the joint distribution of Cj and ^ ■. We have 



/ 5{z - Qdz = 



otherwise 



hence, 



n/2 



j Sn{z)dz = J2 E[cj6jk] = E[ck] =Ck + o{(t] 

Nk ^ = 1 



By a similar argument the second part of the thesis follows. 



□ 
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3 The P-transform 



In order to solve the original moment problem we need to compute 

7l/2 

In order to estimate the expected value we build independent replications of the 
data (pseudosamples) by defining 

^V=^k^^k\ A: = 0,...,n- 1; r=l,...,i? 

where {^'['^''} are i.i.d. zero mean complex Gaussian variables with variance cr'^ in- 
dependent of a/j, V/i. Therefore 

where = cr^ + a'^. For r = 1, . . . , i?, we define the statistics 

n/2 

(r) (r) ■ 

where c^- are the solution of the complex exponentials interpolation problem 

for the data a^^'', /c = 0, . . . , n — 1. As, by Lemma [2], the transformation 

T : {ai'^), k = 0,...,n-l}^ {[cf , ], J = 1, • • • , n/2} 

is one-to-one, Sn,r{z,d''^) are i.i.d. with mean Sn{z,a'^) and finite variance ({z,a'^) 
because {i-'k^} are i.i.d. . Therefore the statistic 

1 ^ . 

has mean Sn{z^d'^) = £"[§^^^(2, 5"^)] and variance ^C,{z^d'^). 
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Let us consider the statistic 

n/2 

Sn(^,0=Ec,<5(z-^,) 

J = l 

where Cj^^j are the solution of the complex exponentials interpolation problem for 
the data a^, A; = 0, . . . , n — 1 and the conditioned statistic 

which are both computable from the observed data a. We have 
Lemma 3 For n and a > fixed and \/z and a, 

lim mr[S^^(z,a2)] = 0. 

proof 

from the conditional variance formula ([23]) we have 
and 

var[{Sl ji{z, a^)] < var[Sn,R{z, a^)] = ^({z, a^). □ 

It follows that \/z the risk of ^(z, a ) as an estimator of S{z) with respect to the 
loss function given by the absolute difference could be smaller than the risk of the 
estimator 8^(2;, cr^) if R and a are suitably chosen, despite of the fact that its bias 
is larger because a > a and Theorem [5] holds. As a matter of fact this possibility is 

33 



always verified provided that a' and R are suitably chosen as proved in the following 



Theorem 6 LetM{z) and Mc{z) be the mean squared error of Sn{z, a'^) an^f ^(z, o"^) 
respectively. In the limit for a ^ 0, it exist a' and R{(j') such that MR > R{(j'), 
M,{z) < M{z) Vz. 

proof 

let Mc{z) = Vc + h1 and M{z) = v + he the decomposition of the mean squared 
errors in the sum of variance plus squared bias. Then M^{z) -b^ = Vc+ {bl - b'^). 
By Lemma El be is equal to the bias of Sn{z,a ) and, by Theorem [5], it is o{a) 
for a ^ 0. Then limo-'^o+(&c ~ = 0- Moreover, by Lemma El limn^ooVc = 0. 
Therefore \/v > 0, 3cr^ and R{(Jy) such that \/a' < a'^, Vc + {bl — b'^) < v and then 
Mc{z) < M{z). □ 

In order to define a discrete transform, we evaluate S^^^(z,a'^)) on a lattice L = 
{(xi, yi),i = 1, . . . , N} such that 

min3?^, > minxj; max3?^, < maxx^ 

j i 3 i 

minQ^o > minyi] maxS'^. < max?/j. 

j i 3 i 

In order to cope with the Dirac distribution appearing in the definition of ^(z, o"^)) 
it is convenient to use an alternative expression given by 



-[ I R n/2 \ 

\.r(^,^') = ( E E[cf |a]log(|^ - Kf|a]|) 



which can be obtained by the former one by remembering that ^Alog(|zp) = 6{z) 
(see e.g. [251 Pg-47]). In this way the problem of discretizing the Dirac 6 is reduced 
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to discretizing the Laplacian operator, which is easier to cope with. We then get a 
random matrix P(a"^) G such that P(/i, k, a"^) = S^^^(x/i + iyk). We call this 

matrix the P-transform of the vector [ag, . . . , SLn^i]. 

4 Estimation procedure 

The P-transform gives a global picture of the measure S{z). However an estimate 
of the unknown parameters {^j, Cj, j = 1, . . . are usually of interest. An auto- 
matic procedure to get such estimates is now described. Let P(cr^) be the P-transform 
computed by using R pseudosamples with variance o"^. The proposed procedure is 
the following (dropping for simplicity the conditioning to a): 

(r) (r) 

• memorize all the Pade' poles and the corresponding residuals , r = 
1, . . . ,R used for computing P(a"^) 

• identify the local maxima of P(5'^) and sort them in increasing order with respect 
to the local maxima values. The local maxima are candidate estimates of {^j, j = 
l,...,p} 

• for each candidate a cluster of (previously memorized) Pade' poles was estimated 
by including all the poles closest to the current candidate until the cluster car- 
dinality equals a predefined percentage (e.g.> 50%) of the number R of pseu- 
dosamples. The rationale is that if the candidate is close to one of the most of 
the pseudosamples should provide a Pade' pole close to it. Notice that spurious 
clusters - i.e. not centered close to some - can be expected [3] 

• all the candidates whose associated cluster does not have the prescribed cardi- 
nality are eliminated. The number p of left candidates is then an estimate of 
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V 

• for each of the p clusters the Pade' poles and the corresponding residuals (previ- 
ously memorized) were then averaged and provided estimates c^, j = 1, . . . ,p 
of the unknown parameters. Hopefully to ij associated to spurious clusters should 
correspond relatively small Cj. 

5 Numerical results 

In this section some experimental evidence of the claims made in the previous 
sections is given. A model with p = 5 components given by 

-0.1-i27r0.3 -0.05-i27r0.28 -0.0001+i27r0.2 -0.0001+i27r0.21 -0.3-i27r0.35" 

c= [6,3,1,1,20], (7 = 0.2, n = 80 

is considered. We notice that SNR = 5 and the frequencies of the 3^^ and 4^^* 
components are closer than the Nyquist frequency (0.21 — 0.20 = 0.01 < 1/n = 
0.0125). Hence a superesolution problem is involved in this case. The quality of 
the approximation of h{z) to the condensed density is first addressed, h{z) is then 
computed along a line which pass through and the closest among the {^h, h ^ j). 
If the model is identifiable h{z) should have a local maximum close to along this 
line. The interquartile range fj of a restriction of h{z) to a neighbor of this maximum 
is then considered as an estimate of the radius of the local support of h{z) assumed 
circular. Then M = 100 independent data sets a^"^^ of length n were generated and 
the Pade' poles <^('"\ m = 1, . . . , M were plotted in figH] where circles of radii f, 
centered on have been represented too. We notice that the circles are reasonable 
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estimates of the Pade' poles clusters which provide an estimate of the support of the 
peaks of the true condensed density corresponding to ^j, j = 1, . . . , j9. We conclude 
that h{z) is a reliable approximation of the condensed density and therefore, with 
the choice of n and a made above, the model is likely to be identifiable. 

We want now to show by means of a small simulation study the quality of the 
estimates of the parameters ^ and c which define the unknown measure S{z). To this 
aim the bias, variance and mean squared error (MSE) of each parameter separately 
will be estimated. M = 500 independent data sets a^™) of length n were generated 
by using the model parameters given above. For m = 1, . . . , M the P-transform 
p("^) was computed based on R = 100 pseudosamples with a'^ = 10~^(7^ on a square 
grid of dimension = 200. The estimation procedure is then applied to each of 
the P^"^\ m = 1, . . . , M and the corresponding estimates cf^\j = 1, . . . of 
the unknown parameters were obtained. As we know the true value if less than p 
local maxima were found in the second step or if p^™^ < p in the fourth step of the 
procedure, the corresponding data set aS"^^ was discarded. 

In Table 1 the bias, variance and MSE of each parameter including p is reported. 
They were computed by choosing among the = 1, . . . the one closest to 

each ^k,k = 1, . . . ,p and the corresponding c^K If more than one is estimated 
by the same ^j"*'' the m— th data set a^'^^ was discarded. In the case considered 65% 
data sets were accepted. Looking at Table 1 we can conclude that the true measure 
can be estimated quite accurately in 65% of cases. 
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When fr^ ' > p we computed also the average residual amplitude 

which represents the contribution to S^^(2;,5-^)) of all the components which give 
rise to spurious clusters. In the case considered its value is ares = 1-165 which 
should be compared with the true amplitudes c. We can conclude that even when 
more components then the true ones are detected their relative importance is very 
low. 

In order to appreciate the advantage of the estimator S^^(2:, a^) with respect to 
8^(2;, (j^), the same M = 100 independent data sets a^^'> of length n generated 
before were considered. The corresponding Pade' poles and weights {ij^\cj^\j = 
1, . . . , n/2) were computed and ordered for each m in decreasing order w.r. to the 
absolute value of the weights. The true {^j,Cj, j = 1, . . . , p) were ordered in the same 
way and the error 

eoH = E(4(-)-^,f +E(cf)-c,f 

i=i j=i 

was computed for m = 1, . . . , M and plotted in fig. 2. Then to each of the M data 
sets cM^^ previously generated R = 100 i.i.d. zero-mean Gaussian samples with 
variance cr'^ = 0.64(7^ were added and Cj^'^\j = 1, . . . , n/2, r = 1, . . . , i?) 

were computed and ordered as before for each m and r. Finally the error 

was computed for m = 1, . . . , M and plotted in fig. 2. We notice that eR{m) <C eo(m) 
for almost all m and it is much less dispersed around its mean. Therefore the 
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estimates of {^j, Cj,j = 1, . . . ,p) obtained by averaging over the R pseudosamples 
are better than those obtained by the original samples. Finally we notice that in 
this simulation we used a variance much larger than the one used to produce the 
results in Table 1. This large value gives the best mean squared error over all the 
five parameters but not necessarily the best reconstruction of each single parameter, 
as we looked for in the previous simulation. 
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P 


bias{p) 


s.d.{p) 


MSE{p) 




5 


0.0500 


1.0000 


1.0025 






bias{ij) 


s.d.ij 


MSEii,) 


J = l 


-0.2796 - 0.8606i 


-0.0006 + 0.0004i 


0.0230 


0.0005 


J = 2 


-0.1782 - 0.9344i 


-0.0005 - 0.0004i 


0.0125 


0.0002 


J = 3 


0.3090 + 0.9510i 


0.0057 - 0.0009i 


0.0171 


0.0003 


J = 4 


0.2487 + 0.9685i 


-0.0005 + 0.0024i 


0.0145 


0.0002 


J = 5 


-0.4354 + 0.5993i 


-0.0054 + 0.0018i 


0.0290 


0.0009 










MSE(cj) 


i = i 


6.0000 


0.1545 


1.7154 


2.9663 


i = 2 


3.0000 


-0.1617 


1.2865 


1.6812 


J = 3 


1.0000 


-0.1037 


0.3295 


0.1193 


J = 4 


1.0000 


-0.0981 


0.3193 


0.1116 


j = 5 


20.0000 


-0.1759 


2.5101 


6.3317 



Table 1 

Statistics of the parameters p, ^j, j = 1, . . . ,p and Cj,j = 1, . . . ,p 
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Fig. 1. Top left: location of Fade' poles for 100 independent realizations of the noise; the circles are the 
estimated support of the condensed density in a neighborhood of ^j; top right:zoom in a neighborhood of 
the 1-st and 2-nd components; bottom left: zoom in a neighborhood of the 3-rd and 4-th components; zoom 
in a neighborhood of the 5-th component (see section 4). 




Fig. 2. MSE of the standard estimator of the parameters i£,j,Cj),j = 1, . . . ,p (dashed); MSE of the averaged 
estimator (solid) 
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